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We present a theory of electronic properties of gated triangular graphene quantum dots with 
zigzag edges as a function of size and carrier density. We focus on electronic correlations, spin and 
geometrical effects using a combination of atomistic tight-binding, Hartree-Fock and configuration 
interaction methods (TB+HF+CI) including long range Coulomb interactions. The single particle 
energy spectrum of triangular dots with zigzag edges exhibits a degenerate shell at the Fermi level 
with a degeneracy N c dge proportional to the edge size. We determine the effect of the electron- 
electron interactions on the ground state, the total spin and the excitation spectrum as a function 
of a shell filling and the degeneracy of the shell using TB+HF+CI for N^dge < 12 and approxi- 
mate CI method for N e d ge > 12. For a half-filled neutral shell we find spin polarized ground state 
for structures up to N = 500 atoms in agreement with previous ab initio and mean-field calcula- 
tions, and in agreement with Lieb's theorem for a Hubbard model on a bipartite lattice. Adding 
a single electron leads to the complete spin depolarization for N e dge < 9- For larger structures, 
the spin depolarization is shown to occur at different filling factors. Away from half-fillings excess 
electrons(holes) are shown to form Wigner-like spin polarized triangular molecules corresponding 
to large gaps in the excitation spectrum. The validity of conclusions is assessed by a comparison 
of results obtained from different levels of approximations. While for the charge neutral system 
all methods give qualitatively similar results, away from the charge neutrality an inclusion of all 
Coulomb scattering terms is necessary to produce results presented here. 



I. INTRODUCTION 

Graphene is an atomically thick layer of carbon atoms 
arranged in a honeycomb lattice. [ll4l0| Due to its unique 
electronic properties and promising potential for appli- 
cations, there is a gro wing re search interest in graphene 
based nanostructures . (lO - 1 2 1 Attempts at fabricating 
graphene nanostructures with well defined shape and 
edge type have been reported starting from the graphene 
layer and usin g to p-down techniques, 13l- 25[ bottom- up 
techniques [26j-|30( starting from carbon based molecules, 
as well as starting from graphane and removing hydrogen 
atoms using AFM tips.|3lH34| 

The work on graphene nanostructures is motivated by 
the expectation that finite size effects significantly mod- 
ify electronic properties of graphene. As a result of size 
quantization, an energy gap opens up, making graphene 
a semiconductor with a gap tunable from THz to UV. 
The energy gap can be tuned by changing not only the 
size but also the shape and the type of edge, allow- 
ing us to control the material's optical properties. [35- 
37[ Two types of edges in graphene are of particular in- 
terest due to their stability: armchair and zigzag. For 
zigzag edges, edge states in the vicinity of the Fermi en- 
ergy appear. This is related to breaking the sublatticc 
symmetry between the two types of atoms in the unit 
cell of the graphene honeycomb lattice. The presence 
of edge states was predicted theoretically [1, 3jJ 38-46] 
and confirmed experimentally. [471 - 1411 These edge states 



form a degenerate band in graphene ribbons 0, 13a - 41 

can collapse to a degenerate shell in graphene quantum investigated. [4 



dots. [35L l42T - |46l . [Hol. |51 1 It was previously shown that the 
degeneracy is equal to the difference between the number 
of atoms corresponding to two sublattices in the bipartite 
lattice. 0, [43], 0, [H3| In particular, the geometry that 
maximizes the imbalance between the two sublattices is 
a zigzag edge triangle where the degeneracy of the zero- 
energy shell is proportional to the number of atoms on 
the one edge. [46} This presents a unique opportunity to 
design a quantum system with a macroscopic degener- 
acy, analogously to the two-dimensional electron gas in a 
strong magnetic field. 

While fabricating and measuring triangul ar grap hene 
quantum dots with well defined edges |2C 



remains a challenge, the theory of triangular graphene 
quantum dots (TGQD) with z igzag e dges was developed 
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by several groups. [29|, |35j, |3JJ, I42M4 
ticular, the macroscopically degenerate zero-energy band 
and the corresponding wavefunctions were explicitely 
constructed. [46j For a half- filled shell, TGQDs were stud- 
ied by Ezawa using the Heisenberg Hamiltonian, [42I ] by 
Fernandez-Rossier and Palacios [43J using the mean-field 
Hubbard model, by Wang, Meng, and Kaxiras using 



density functional theory (DFT); and Giiclii et al. [51 1 
using exact diagonalization techniques. It was shown 
that the ground state is fully spin polarized, with a fi- 
nite magnetic moment proportional to the shell degen- 
eracy. This finding is in agreement with Lieb's theorem 
on magnetism of the Hubbard model for bipartite lattice 
systems. 65 

The effect of defects and disorder was also 
571 . l58l ] In particular, Voznyy et 



al. 
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[58j have shown by using ab initio methods that 
hydrogen-passivation stabilizes zigzag edges in TGQD 
over the pentagon- heptagon reconstructed edges. [58j| 
It was also proved that the zero-energy shell survives 
when TGQD is deformed to trapezoidal shape, (iij 
Ezawa studied the stability of magnetization against 
disorder. He considered three types of randomness: in 
a hopping integral, a site energy and lattice defects. [57j 
He proved that the magnetism is still governed by Lieb's 
theorem but the number of degenerate states changed 
by the number of lattice defects. Some of us have 
shown in Ref. 51 by use of methods beyond mean- field 
approximations, that the magnetization is unstable with 
respect to additional charge, leading to a complete spin 
depolarization. The spin depolarization was shown to 
significantly influence transport properties, blocking 
current through the graphene quantum dot due to the 
spin blockade. [51| It was also shown that by changing 
the population of the degenerate shell using a gate, 
one can simultaneously control magnetic and optical 
properties, determined by strong electron-electron and 
excitonic interactions. 
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In this work we use improved configuration-interaction 
(CI) tools to extend our previous results [5l| regarding 
the role of electron-electron interactions, magnetism and 
correlations in TGQDs to larger structures. We investi- 
gate the electronic properties as a function of size and fill- 
ing factor of the degenerate shell by using a combination 
of tight-binding (TB), Hartree-Fock (HF) and configura- 
tion interaction methods (TB+HF+CI). Our many-body 
Hamiltonian includes, in addition to the on-site inter- 
action term, all scattering and exchange terms within 
next-nearest neighbors, and all direct interaction terms 
in the two-body Coulomb matrix elements. Using full 
CI combined with the TB+HF method we demonstrate 
that the ground state for the charge neutral system has 
maximally polarized edge states for structures consisting 
of up to 200 atoms with the number of degenerate edge 
states N e d ge < 9. By analyzing a spin-flip excitation 
spectrum of the spin-polarized ground state, we verify 
the spin-polarized ground state for up to 500 atoms or 
N e dge = 20. These results for a system with long ranged 
Coulomb interaction appear to be consistent with Lieb's 
theorem for the Hubbard model. Using TB+HF+full CI 
method for TGQD charged with an additional electron 
and a size of up to N — 200 atoms it is shown that a 
complete spin depolarization predicted earlier by some 
of us [5l| exists only up to a critical size. The critical 
size is established by studying the stability of a charged 
spin-polarized shell to spin-flip excitations. It is shown 
that for sizes up to the critical size the spin wave and mi- 
nority spin electron form a bound state, a trion, signaling 
the tendency to the depolarization. For sizes exceeding 
the critical size the spin waves are unbound and the spin- 
polarized state is the ground state up to the sizes studied 
(N sa 500 atoms). For TGQD structures above the crit- 



ical size, depolarization effects away from the half-filling 
are observed. Results of TB+HF+CI calculations allow 
us to extract the excitation gap as a function of a shell 
filling. It is found that the largest gaps correspond to 
the half-filled spin-polarized shell and special filling frac- 
tions. At these special filling fractions, we predict a for- 
mation of Wigner-like spin polarized molecules, related 
to long range Coulomb interactions and a triangular ge- 
ometry of graphene quantum dot. Finally, we compare 
results obtained at different levels of approximations. We 
show that, for the charge neutral system, the Hubbard, 
extended Hubbard, and fully interacting models are in 
good qualitative agreement. On the other hand, away 
from the half-filling, only a fully interacting model is able 
to correctly capture the effect of correlations. 

The paper is organized as follows. In Sec. II, we de- 
scribe our model. Section III contains analysis of the 
ground state spin and correlations as a function of size 
and filling factor of the degenerate shell. In Sec. IV, 
we compare results obtained within different levels of ap- 
proximations. In Sec. V, we summarize our results. 



II. MODEL OF A GRAPHENE TRIANGULAR 
QUANTUM DOT 

Graphene is a two-dimensional honeycomb crystal 
formed by carbon atoms on two interpenetrating hexag- 
onal sublattices. The unit cell thus contains two carbon 
atoms. The distance between nearest-neighbor atoms 
is a = 1.42 A. By using vectors R = nai + m&2 
with n, m integers and primitive unit vectors defined as 
a-i.2 = a/2(±V3, 3), one can obtain the positions of all 
the atoms in the structure. By cutting the graphene lat- 
tice in three zigzag directions, an equilateral triangle can 
be obtained, as shown in Fig. Q] Such a system has a bro- 
ken sublattice symmetry with two properties: (i) all edge 
atoms (with only two bonds) belong to the same sublat- 
tice, (ii) the difference between the number of atoms be- 
longing to each sublattice is proportional to the number 
of atoms on one of the three edges. 

Each carbon atom has four valence electrons. Three 
of them, on s, p x , and p y orbitals, form sp 2 bonds 
with the three nearest in-plane neighbors. They are 
strongly bound and responsible for mechanical proper- 
ties of graphene. The remaining fourth valence electron 
on each carbon atom p z orbital, perpendicular to the 
plane of graphene, is weakly bound and determines elec- 
tronic properties of the system. Single-particle properties 
of graphene can be described by using one orbital tight- 
binding (TB) Hamiltonian. [Hij We have previously shown 
that, within the TB model in the nearest- neighbors ap- 
proximation, TGQDs with zigzag edges exhibit an en- 
ergy gap, with a degenerate shell at the Fermi (zero) en- 
ergy, with a degeneracy proportional to the length of an 
edge. 4G] An example of TB energy levels for a structure 
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FIG. 1: (Color online) Triangular graphene quantum dot with 
zigzag edges. There are eight edge atoms (with two bonds) 
on one edge. Red (light gray) and blue (dark gray) colors dis- 
tinguish between two sublattices in the honeycomb graphene 
lattice. Structure consists of a total of N = 97 atoms 



consisting of 97 atoms with N e d ge = 7 degenerate states 
is shown in Fig. Our goal is to study the role of 

electron-electron interactions for electrons occupying this 
degenerate shell. Solving the full many-body problem 
even for such a small structure with 97 atoms is not pos- 
sible at present. However, due to the energy gap separat- 
ing the valence band and degenerate states, the valence 
electrons that do not occupy the degenerate band can be 
treated in a mean-field approximation. The remaining 
electrons occupying the degenerate shell must, however, 
be treated using a configuration-interaction method (CI). 
Therefore, we use a TB+HF+CI approach that allows us 
to treat the electronic correlations for electrons in the de- 
generate shell and their interaction with valence electrons 
at the mean-field level. 

We start from the full many-body Hamiltonian for in- 
teracting electrons on the p z orbitals of graphene. It can 
be written as 



> 

LU 




44 46 48 50 

eigenstate index 



52 



FIG. 2: (a) Single-particle nearest-neighbor tight-binding 
(TB) energy levels. The zero-energy shell on the Fermi level 
is perfectly degenerate, (b) Positively charged system with an 
empty degenerate band after self-consistent Hartee-Fock (HF) 
mean-field calculations described by a single Slater determi- 
nant (TB+HF model), (c) Occupation of empty degenerate 
HF quasi-orbitals by electrons. The inset pictures schemat- 
ically show the excess charge corresponding to each of the 
three model systems. The ground state and the total spin of 
the system of interacting electrons can be calculated by us- 
ing the configuration interaction (CI) method, described in 
Sec. II (TB+HF+CI). The charge neutrality corresponds to 
a half- filled degenerate band (not shown). 



A. Mean-Field approximation for infinite graphene 
sheet 



Let us first write the Hamiltonian for graphene, given 
by Eq. ([1]), in the mean- field Hartree-Fock (HF) approx- 
imation as 



H = E T il a c\ a Cl a + - E (^\ V \ kl )4a C ]a' C k<y'Cla, (1) 
go 1 



where the operator c i(J creates an electron on the i-th p z 
orbital with spin cr, Tu a is a hopping integral that de- 
scribes the probability of scattering of electron on the 
Z-th p z orbital fa to the i-th p z orbital fa. The sec- 
ond term describes two-body Coulomb interactions be- 
tween p z electrons. Note that at this stage, the unknown 
hopping terms Tua do not include the effect of electron- 
electron interactions of p z electrons. 



h°mf = E^V + E E p° ka >mv\ki) 

- {ij\V\lk)5^,)clc la = E Ui*cic la . (2) 

This is effectively a one-body TB Hamiltonian for a 
graphene layer [66( with density matrix elements p° ka , — 

( c ] a ' c ka') calculated with respect to the fully occupied 
valence band. The values of p° k(J , are evaluated in Ap- 
pendix A and their role becomes clear in the next subsec- 
tion, tua are experimentally estimated hopping integrals. 
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B. Mean-Field approximation for graphene quantum 

dots 



We now derive a mean-field Hamiltonian for electrons 
in graphene quantum dots (GQD). First, we apply mean- 
field approximation to the Hamiltonian given by Eq. ([1]) 
for electrons in GQD, with a result written as 



H 



GQD 
MF 



\ nia^cia + 2J ^2 P^ ka ' i(^\ V \ kl ) densities shown in Ref. El 



integral tua- On the other hand, close to the edges, 
a density matrix for the GQD differs in comparison to 
its graphene counterpart. After diagonalizing the HF 
Hamiltonian given by Eq. Q one obtains eigenvalues 
and eigenvectors that involve the geometrical properties 
of the system, shown in Fig. [U[b) • A slight removal of the 
degeneracy of middle edge states and three corner states 
with a bit higher energies are observed, with electronic 



{ij\V\lk)5 a ^)c\ a ci a , 



(3) 



with density matrix pjka- 1 for GQD. By combining Eq. 
© and Eq. (0]) we get 



ttGQD jjGQD n-o . tjo 

n MF — MF ~ a MF T n MF 



- {ij\V\lk)8 a y)clc lB 



(4) 



Here the subtracted component in the second term cor- 
responds to mean-field interactions included in effective 
Uia hopping integrals, described by the graphene density 
matrix p° ka ,- For the TGQDs, the density matrix ele- 
ments Pjka' are calculated with respect to the many-body 
ground state of N re f = N sit e — N e d ge electrons, where 
N S ite is the number of atoms. Since the valence band 
and the degenerate shell are separated by an energy gap, 
the closed shell system of N re f interacting electrons is 
expected to be well described in a mean-field approxima- 
tion, using a single Slater determinant. This corresponds 
to a charged system with N e d ge positive charges, as 
schematically shown in Fig. EJb). The Hamiltonian given 
by Eq. ^ has to be solved self-consistently to obtain 
Hartree-Fock quasi particle orbitals. In numerical cal- 
culations, in addition to the on-site interaction term, all 
scattering and exchange terms within next-nearest neigh- 
bors and all direct interaction terms are included in the 
two-body Coulomb matrix elements (ij\V\kl) computed 
using Slater p z orbitals. 67] The few largest Coulomb ma- 
trix elements are given in Ref. |68|. The value of the effec- 



tive dielectric constant k depends on the substrate and 
is set to k — 6 in what follows. [H^| A method of calcu- 
lating values of p° ka , for graphene is shown in Appendix 
A. Values of hopping integrals tu a are taken from the 
experimental data 70] or ab initio calculation. [69j We 
use t = —2.5 eV for nearest-neighbors and t' = —0.1 eV 
for next-nearest-neighbors [rij j hopping matrix elements. 
The HF results were also compared with the results of ab 
initio calculations. 
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We now discuss mean-field results for the charge neu- 
tral system. In the vicinity of the center of a sufficiently 
large dot a charge distribution around a given site is iden- 
tical to that of an infinite system. The density matrices 
for graphene p° ka , and for GQD pjka' are equal. A sec- 



C. Configuration Interaction method 

After the self-consistent procedure we get new orbitals 
for quasi-particles with a fully occupied valence band and 
a completely empty degenerate shell. We start filling 
these degenerate states by adding extra electrons one by 
one, schematically shown in Fig. EJc). Next, we solve 
the many-body Hamiltonian corresponding to the added 
electrons, given by 

Hmb — ^ £ s al a a sa 

s,<7 

+ \ (sp I V I df)al a al <J ,a da 'a fa , (5) 

where the first term describes the energies of the Hartree- 
Fock orbitals and the second term describes an interac- 
tion between quasi particles occupying degenerate HF 
states denoted by s,p,d,f indices. The two-body quasi 
particle scattering matrix elements (sp \ V \ df) are cal- 
culated from the two-body localized on-site Coulomb ma- 
trix elements (ij \V\ kl). 

In our calculations, we neglect scattering from/to the 
states from a fully occupied valence band. Moreover, 
because of the large energy gap between the shell and 
the conduction band, we can neglect scatterings to the 
higher energy states. A validity of these approximations 
is assessed in Ref. 68]. These approximations allow us 
to treat the degenerate shell as an independent system 
that significantly reduces the dimension of the Hilbert 
space. The basis is constructed from vectors correspond- 
ing to all possible many-body configurations of electrons 
distributed within the degenerate shell. For a given num- 
ber of electrons N e i , the Hamiltonian given by Eq. is 
diagonalized in each subspace with total S z . 

D. Effect of gate charge 

In our model, we start from the system with an empty 
shell that corresponds to the charged system. As in our 
previous work, Ref. |5jJ, electrons from the shell are 
transferred to the metallic gate. The Hamiltonian for 



ond term in Eq. Q vanishes, leaving only a hopping N re f electrons in the presence of a gate in the mean-field 
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Hartree-Fock approximation was written as 
H MF = tu^lcia + ^ k -' - P°jk^)i(ij\V\kl) 

- (ij\V\lk)5 !7 ,a')c l la ci a + J2vf l (q lnd )cl (y c la , (6) 

i,a 

with an electrostatic potential related to N edge elec- 
trons in a gate defined as 



N eit , 

vfiiqind) = 



-qind/N s ite 



<7) 



~i k^J {xi - Xj) 2 + (y t - yj) 2 + d 2 ga1 

with qi n( i = —N e d ge charge smeared out at positions 
(xi,yi) at a distance d ga te from the quantum dot. Next, 
we derived the many-body Hamiltonian with an inclusion 
of the effect of gate, written as 

a. a 1 

+ J2 (p\v 9 (N add )\q)al a a qa + 2 £<PV(^«w)b'>, (8) 

where the indices without the prime sign (p, q, r, s) run 
over N edge degenerate states, while the index with the 
prime sign p' runs over 7V re //2 valence states (below the 
degenerate shell). A third term in Eq. (J8)) corresponds to 
scattering from state q to state p in a degenerate shell as a 
result of interactions with electrons in a gate. The fourth 
term is a constant and just shifts the entire spectrum by 
a constant energy. 



III. MAGNETISM AND CORRELATION 
EFFECTS 



Electronic properties as a function of the filling 
factor 



We, first, concentrate on a TGQD consisting of N = 97 
atoms, which is the largest system previously studied in 
our earlier work in Ref. [5l|. It has N edge — 7 zero- 
energy degenerate states obtained from TB calculations, 
shown in Fig. HJa). After self-consistent HF calculations 
neglecting the gate charge (the effect of the gate will be 
discussed later), we obtain new quasi particle orbitals, 
shown in Fig. [2jb) . The degeneracy is slightly removed. 
We fill these degenerate levels by additional electrons and 
calculate two-body scattering matrix elements. For a 
given number of quasi particles, the many-body Hamilto- 
nian, Eq. ([5]), is diagonalized in a basis of configurations 
of electrons distributed within the shell, as explained in 
Sec. II. In Fig. El we analyze the dependence of the 
low energy spectra on the total spin S for [Fig. E[a)] 
the charge neutral system, N e i = 7 electrons, and [Fig. 
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FIG. 3: The low-energy spectra for the different total spin 
S for (a) N e i — 7 electrons and (b) N e i = 8 electrons. For 
N e i = 7 electrons the ground state corresponding to S = 3.5, 
indicated by a circle, is well separated from excited states 
with different total spin 5*. For N e i — 8 electrons the ground 
state corresponding to S = 0, indicated by a circle, is almost 
degenerate with excited states with different total spin 5. 



[3fb)] one added electron, i.e., 7V e ; = 8 electrons. We see 
that for the charge neutral TGQD with N e i = 7 elec- 
trons the ground state of the system is maximally spin 
polarized, with S — 3.5, indicated by a circle. There is 
only one possible configuration of all electrons with par- 
allel spins that corresponds to exactly one electron per 
one degenerate state. The energy of this configuration 
is well separated from other states with lower total spin 
S, which require at least one flipped spin among seven 
initially spin-polarized electrons. An addition of one ex- 
tra electron to the system with N e i — 7 spin polarized 
electrons induces correlations as seen in Fig. [3jb) , where 
the cost of flipping one spin is very small. Moreover, for 
N e i = 8, the ground state is completely depolarized with 
5 = 0, indicated by a circle, but this ground state is 
almost degenerate with states corresponding to the dif- 
ferent total spin. 

The calculated many-body energy levels, including all 
spin states for different numbers of electrons (shell fill- 
ing), are shown in Fig. HI For each electron number, 
N e i , energies are measured from the ground-state energy 
and scaled by the energy gap of the half-filled shell, cor- 
responding to N e i = 7 electrons in this case. The solid 
line shows the evolution of the energy gap as a func- 
tion of shell filling. The energy gaps for a neutral sys- 
tem , N e i — 7 , as well as for N e i = 7 — 3 = 4 and 
N e i = 7 + 3 = 10 are found to be significantly larger in 
comparison to the energy gaps for other electron num- 
bers. In addition, close to the half-filled degenerate shell, 
the reduction of the energy gap is accompanied by an 
increase of low energy density of states. This is a signa- 
ture of correlation effects, showing that they can play an 
important role at different filling factors. 

We now extract the total spin and energy gap for each 
electron number. Figures [5ja) and (b) show the phase 
diagram, the total spin S and an excitation gap as a 
function of the number of electrons occupying the degen- 



FIG. 6: (Color online) The spin densities of the ground state 

0.0- for (a) N e i = 4 electrons and (b) N e i = 10 electrons that 

correspond to subtracting/adding three electrons from/to the 
1 2 3 4 5 6 7 8 91011121314 charge neutral system. The radius of circles is proportional to 

^ a value of spin density on a given atom. A long range Coulomb 

el interaction repels (a) holes and (b) electrons to three corners, 

forming a spin-polarized Wigner-like molecule. 



FIG. 4: The low-energy spectra of the many-body states as a 
function of the number of electrons occupying the degenerate 
shell for the system with N e d g e = 7 degenerate states. The 
energies are renormalized by the energy gap corresponding 
to the half-filled shell, N e i — 7 electrons. A large density of 
states around N e d ge + 1 electrons is a signature of the correla- 
tion effects. The solid line shows the evolution of the energy 
gap as a function of shell filling. 
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FIG. 5: (Color online) (a) The total spin as a function of 
number of electrons occupying the degenerate shell and (b) 
corresponding the energy excitation gaps, with and without a 
gate, red (light gray) and black (dark gray) lines, respectively. 
Due to a presence of correlation effects for some fillings, the 
magnitude of the energy gap is significantly reduced. 



erate shell. The system reveals maximal spin polarization 
for almost all fillings, with exceptions for N e i =8,9 elec- 
trons. However, the energy gaps are found to strongly 
oscillate as a function of shell filling as a result of a com- 
bined effect of correlations and system's geometry. We 
observe a competition between fully spin polarized sys- 
tem that maximizes exchange energy and fully unpolar- 
ized system that maximizes the correlation energy. Only 
close to the charge neutrality, for N e i = 8 and N e i = 9 
electrons, are the correlations sufficiently strong to over- 
come the large cost of the exchange energy related to 
flipping spin. The excitation gap is significantly reduced 
and exhibits large density of states at low energies, as 
shown in Fig. [3] 

Away from half-filling, we observe larger excitation 
gaps for N e i = 4 and N e i = 10 electrons. These fill- 
ings correspond to subtracting/adding three electrons 
from/to the charge- neutral system with N e i = 7 elec- 
trons. In Fig. IS] we show the corresponding spin densi- 
ties. Here, long range interactions dominate the physics 
and three spin polarized [Fig. HJa)] holes (N e i =7 — 3 
electrons) and [Fig. [Bfb)] electrons (N e j = 7+3 electrons) 
maximize their relative distance by occupying three con- 
secutive corners. Electron spin density is localized in 
each corner while holes correspond to missing spin den- 
sity localized in each corner. We also note that this is not 
observed for N e i = 3 electrons filling the degenerate shell 
(not shown here). The energies of HF orbitals of corner 
states correspond to three higher energy levels [see Fig. 
H^lc)]. with electronic densities shown in Ref. [5l|. Thus, 
N e i = 3 electrons occupy lower-energy degenerate levels 
corresponding to sides instead of corners. On the other 
hand, when N e i = 7 electrons arc added to the shell, self- 
energies of extra electrons renormalize the energies of HF 
orbitals. The degenerate shell is again almost perfectly 
flat similarly to levels obtained within the TB model. A 
kinetic energy does not play a role allowing a formation 
of a spin-polarized Wigner-like molecule, resulting from 
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Electronic properties as a function of the size 




d , [a] 

gate L J 



FIG. 7: The energy gaps around the charge neutrality for a 
system with N e d ge = 7 degenerate states as a function of a 
gate distance. The energy gap for the charge-neutral system, 
N e i = 7, changes by less than 1%. For the charged system, 
N e i — 8, we observe changes in the energy gap for a gate 
distance in a range 5a < d ga te < 10a but still not affecting 
the spin depolarization. 



a long-range interactions and a triangular geometry. We 
note that Wigner molecules were previously discussed in 
circular graphene quantum dots with zigzag edges de- 
scribed in the effective mass approximation. 72|, [73| The 
rotational symmetry of quantum dot allowed for the con- 
struction of an approximate correlated ground state cor- 
responding to either a Wigner-crystal or Laughlin-likc 
state. [72[ Later, a variational rotating-electron-molecule 
(VREM) wave function was used.[73j Unfortunately, due 
to a lack of an analytical form of a correlated wave func- 
tion with a triangular symmetry, it is not possible to do 
it here. 

Figure [5] also shows the effect of the presence of a gate 
at a distance d ga t e = 20a, where a = 1.42 is a nearest- 
neighbor's inter-atomic distance. Clearly, the effect of 
a gate is very weak, just slightly changing energy gaps. 
In Fig. [71 energy gaps as a function of a gate distance 
for the charge-neutral N e i — 7 and charged N e i = 8 sys- 
tem for our tested system with N e d ge = 7 degenerate 
states are shown. There are no effects for a gate dis- 
tance dgate > 20a. When a gate distance is comparable 
to graphene-substrate separation, d ga t e ~ 5a, the energy 
gap for N e i = 7 increases while the energy gap for N e i = 8 
decreases. The drop for N e i — 8 is not sufficiently strong 
to change an observed effect of the spin depolarization. 
According to the above analysis, we next present results 
for a Hamiltonian with a gate at infinity. 



In a previous section, we have analyzed in detail the 
electronic properties of a particular TGQD with N = 97 
atoms as a function of the filling factor v — N e i/N e d ge , 
i.e., the number of electrons per number of degenerate 
levels. In this section we address the important question 
of whether one can predict the electronic properties of a 
TGQD as a function of size. 

Figure [8] shows spin phase diagrams for triangles with 
odd number of degenerate edge states N e d ge and increas- 
ing size. Clearly, the total spin depends on the filling fac- 
tor and size of the triangle. However, all charge-neutral 
systems at v = 1 are always maximally spin polarized 
and a complete depolarization occurs for N e d ge < 9 for 
structures with one extra electron added (such depolar- 
ization also occurs for even N e d ge , not shown). Similar 
results for small size triangles were obtained in our previ- 
ous work.[5l( However, at N e d ge = 11 we do not observe 
depolarization for N e d ge + 1 electrons but for N e d ge + 3, 
where a formation of Wigner-like molecule for a triangle 
with N e dge = 7 was observed. We will come back to this 
problem later. We now focus on the properties close to 
the charge neutrality. 

For the charge-neutral case, the ground state corre- 
sponds to only one configuration \GS) — \\ i aj JO) with 
maximum total S z and occupation of all degenerate shell 
levels i by electrons with parallel spin. Here |0) is the HF 
ground state of all valence electrons. Let us consider the 
stability of the spin polarized state to single spin flips. 
We construct spin-flip excitations \kl) = a\ 1 .ai^\GS) 
from the spin-polarized degenerate shell. The spin-up 
electron interacts with a spin-down "hole" in a spin- 
polarized state and forms a collective excitation, an ex- 
citon. An exciton spectrum is obtained by building an 
exciton Hamiltonian in the space of electron- hole pair ex- 
citations and diagonalizing it numerically, as was done, 
e.g., for quantum dots. 74| If the energy of the spin flip 
excitation turns out to be negative in comparison with 
the spin-polarized ground state, the exciton is bound and 
the spin-polarized state is unstable. The binding energy 
of a spin-flip exciton is a difference between the energy 
of the lowest state with S — S" lax — 1 and the energy 
of the spin-polarized ground state with S = S™ ax . An 
advantage of this approach is the ability to test the sta- 
bility of the spin polarized ground state for much larger 
TGQD sizes. 

Figure[9]shows the exciton binding energy as a function 
of the size of TGQD, labeled by a number of the degen- 
erate states N e d ge - The largest system, with N e d ge = 20, 
corresponds to a structure consisting of N = 526 atoms. 
The exciton binding energies are always positive, i.e., 
the exciton does not form a bound state, confirming a 
stable magnetization of the charge neutral system. The 
observed ferromagnetic order was also found by other 
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FIG. 8: Spin phase diagrams as a function of filling fac- 
tor v — N e i/N e dge for different size triangles characterized by 
the number of the degenerate edge states Nedge- Half- filled 
shell v = 1 is always maximally spin polarized. The com- 
plete spin depolarization occurs for one added electron to the 
charge-neutral system for Nedge < 9. For N e d ge = 11 the 
depolarization effect moves to a different filling. 



groups based on calculations for small systems with dif- 
ferent levels of approximations. [42, 43, 50, El| The above 
results confirm predictions based on Lieb's theorem for a 
Hubbard model on bipartite lattice relating total spin to 
the broken sublattice symmetry. [65| Unlike in Lieb's the- 
orem, in our calculations many-body interacting Hamil- 
tonian contains direct long-range, exchange, and scatter- 
ing terms. Moreover, we include next-nearest-neighbor 
hopping integral in HF self-consistent calculations that 
slightly violates bipartite lattice property of the system, 
one of cornerstones of Lieb's arguments. [65J Nevertheless, 
the main result of the spin-polarized ground state for the 
charge neutral TGQD seems to be consistent with pre- 
dictions of Lieb's theorem and, hence, applicable to much 
larger systems. 

Having established the spin polarization of the charge- 
neutral TGQD we now discuss the spin of charged 
TGQD. We start with a spin-polarized ground state \GS) 
of a charge-neutral TGQD with all electron spins down 
and add to it a minority spin electron in any of the 
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FIG. 9: Size-dependent analysis based on exciton and trion 
binding energies. For the charge-neutral system, it is energet- 
ically unfavorable to form an exciton, which is characterized 
by a positive binding energy. Observed dependence confirms 
Lieb's theorem regarding the magnetization of the bipartite 
lattice systems. The formation of a trion is desirable for small 
size systems. The phase transition occurs close to N e dge = 8, 
indicated by an arrow. The complete depolarization effects 
close to the charge neutrality observed previously in TGQD 
with N e dge < 9 for Nedge + 1 electrons in Fig. [8] is predicted 
to not appear for larger systems. 



spin of these states is S™ ax — 1/2. We next study sta- 
bility of such states with one minority spin-up electron 
to spin-flip excitations by forming three particle states 
\lki) = al t a k>i al t \GS) with total spin S™ ax -1/2-1. 
Here there are two spin-up electrons and one hole with 
spin-down in the spin-polarized ground state. The inter- 
action between the two electrons and a hole leads to the 
formation of trion states. We form a Hamiltonian matrix 
in the space of three particle configurations and diago- 
nalize it to obtain trion states. If the energy of the lowest 
trion state with S™ ax — 1/2 — 1 is lower than the energy 
of any of the charged TGQD states \i) with S™ ax - 1/2, 
the minority spin electron forms a bound state with the 
spin-flip exciton, a trion, and the spin-polarized state of 
a charged TGQD is unstable. The trion binding energy, 
shown in Fig. [51 is found to be negative for small sys- 
tems with N ec ige < 8 and positive for all larger systems 
studied here. The binding of the trion, i.e., the nega- 
tive binding energy, is consistent with the complete spin 
depolarization obtained using TB+HF+CI method for 
TGQD with N e d g e < 9 but not observed for N e d g e = 11 
(and not observed for N e dge = 10, not shown here), as 
shown in Fig. [8] For small systems, a minority spin-up 
electron triggers spin-flip excitations, which leads to the 
spin depolarization. With increasing size, the effect of 
the correlations close to the charge neutrality vanishes. 
At a critical size, around N e d ge — 8, indicated by an ar- 
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FIG. 10: The low-energy spectra of the many-body states as a 
function of the number of electrons occupying the degenerate 
shell for the triangle with N e dge — 11 degenerate states. The 
energies are renormalized by the energy gap corresponding to 
the half-filled shell, N e i = 11 electrons. The large density 
of states related to the correlation effects observed in Fig. [4] 
around N e d g e + 1 electrons shifts to a different filling around 
N e d g e + 3 electrons. 



row in Fig. [§1 a quantum phase transition occurs from 
minimum to maximum total spin. 

However, the spin depolarization does not vanish but 
moves to different filling factors. In Fig. [8] we observe 
that the minimum spin state for the largest structure 
computed by the TB+HF+CI method with N edge = 11 
occurs for TGQD charged with additional three elec- 
trons. We recall that for TGQD with N e dge = 7 charged 
with three additional electrons a formation of a Wigner- 
like spin polarized molecule was observed, shown in Fig. 
IH1 In the following, the differences in the behavior of 
these two systems, N e dge — 7 and 11, will be explained 
based on the analysis of the many-body spectrum of the 
Nedge = 11 system. 

Figure [TU] shows the many-body energy spectra for dif- 
ferent numbers of electrons for N e d g e = H TGQD to 
be compared with Fig. @] for the N e d ge = 7 structure. 
Energies are renormalized by the energy gap of a half- 
filled shell, N e i = 11 electrons in this case. In contrast to 
the N e dg e = 7 structure, energy levels corresponding to 
N e i = Nedge + 1 electrons are sparse, whereas increased 
low-energy densities of states appear for N e j = N e dge + 2 
and N e i — N e dge + 3 electrons. In this structure, elec- 
trons are not as strongly confined as for smaller sys- 
tems. Therefore, for N e i — N e dge + 3 electrons, geo- 
metrical effects that lead to the formation of a Wigner- 
like molecule become less important. Here, correlations 
dominate, which results in a large low-energy density of 
states. 




FIG. 11: Hartree-Fock energy levels corresponding to the de- 
generate shell for calculations with (a) only the on-site term 
U (Hubbard model), (b) the on-site term U + direct long 
range interaction (extended Hubbard model), and (c) all in- 
teractions. A separation of three corner states with higher 
energies is related to direct long range Coulomb interaction 
terms. 



IV. DIFFERENT LEVELS OF 
APPROXIMATIONS ANALYSIS 

In this section, we study the role of different interaction 
terms included in our calculations. The computational 
procedure is identical to that described in Sec. II. We 
start from the TB model but in self-consistent HF and 
CI calculations we include only specific Coulomb matrix 
elements. We compare results obtained with Hubbard 
model with only the on-site term, the extended Hub- 
bard model with on-site plus long range Coulomb inter- 
actions, and a model with all direct and exchange terms 
calculated for up to next-nearest neighbors using Slater 
orbitals, and all longer range direct Coulomb interaction 
terms approximated as (ij\V\ji) = l/(K,\ri — rj\), written 
in atomic units, 1 a.u.= 27.211 eV, where r< and rj are 
positions of i-th and j-th sites, respectively. 

The comparison of HF energy levels for the structure 
with Nedge = 7 is shown in Fig. QTJ The on-site U- 
tcrm slightly removes degeneracy of the perfectly fiat 
shell [FigfTTTa)] and unveils the double valley degeneracy. 
On the other hand, the direct long-range Coulomb inter- 
action separates three corner states from the rest with a 
higher energy [Fig fTTT b)]. forcing the lifting of one of the 
doubly degenerate subshells. Finally, the inclusion of ex- 
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FIG. 12: Spin phase diagrams obtained by use of the CI 
method with (a) only the on-site term U (Hubbard model), 
(b) the on-site term U + direct long range interaction (ex- 
tended Hubbard model), and (c) all interactions. The fer- 
romagnetic order for the charge-neutral system is properly 
predicted by all three methods. Correlations leading to the 
complete depolarization for N e i = Nedge + 1 electrons and 
N e i = N e d a e + 2 electrons are observed only within a full 
interacting Hamiltonian. 



change and scattering terms causes stronger removal of 
the degeneracy and changes the order of the four lower- 
lying states. However, the form of the HF orbitals is not 
affected significantly (not shown here). 

In Fig. Q2] we study the influence of different inter- 
acting terms on CI results. The phase diagrams ob- 
tained within (a) the Hubbard model and (b) the ex- 
tended Hubbard model show that all electronic phases 
are almost always fully spin polarized. The ferromag- 
netic order for the charge-neutral system is properly pre- 
dicted. For TGQD charged with electrons, only inclu- 
sion of all Coulomb matrix elements correctly predicts 
the effect of the correlations leading to the complete de- 
polarization for N e i — 8 and 9. We note that the de- 
polarizations at other filling factors are also observed in 
Hubbard (at N e i = 2)) and extended Hubbard calcula- 
tion (at N e i = 11)) results. 



A more detailed analysis can be done by looking at 
the energy excitation gaps, which are shown in Fig. 1131 
For the charge-neutral system, all three methods give 
comparable excitation gaps, in agreement with previous 



results. [42], |43|, l5fj, l5l| In the Hubbard model, the en- 



ergy gap of the doped system is reduced compared to the 
charge neutrality but without affecting magnetic proper- 
ties. The inclusion of a direct long-range interaction in 
Fig. fTBTb) induces oscillations of the energy gap. For 
N e i = N e dge + 1 electrons the energy gap is significantly 
reduced but the effect is not sufficiently strong to de- 
polarize the system. Further away from half-filling, a 
large energy gap for models with long-range interactions 
for N e i — Nedge + 3 appears, corresponding to the for- 
mation of a Wigner-like molecule of three spin-polarized 
electrons in three different corners. The inclusion of ex- 
change and scattering terms slightly reduces the gap but 
without changing a main effect of Wigner-like molecule 
formation. 



V. CONCLUSIONS AND REMARKS 

We have investigated magnetism, correlations, and ge- 
ometrical effects in TGQDs by use of the TB+HF+CI 
method. Our many-body Hamiltonian includes all direct 
long-range terms and exchange and scattering terms up 
to next-nearest neighbors. We have performed analysis 
as a function of the filling factor of the degenerate band 
of edge states for different sizes. Through a full analysis 
of the many-body energy spectrum of structures consist- 
ing of up to 200 atoms, we confirmed the existence of 
the spin polarized ground state in agreement with Lieb's 
theorem. By studying spin exciton binding energies, we 
also predicted stable magnetization for structures with 
more than 500 atoms. The complete spin depolarization 
was observed for one electron added to the charge neutral 
TGQD up to a critical size. Above a critical size the max- 
imally spin-polarized charged TGQD was predicted us- 
ing trion binding energy analysis. We have shown that in 
small systems, three electrons/holes added to the charge 
neutrality form the spin-polarized Wigner-like molecule. 
We relate this fact to geometrical effects and direct long- 
range interaction terms. For larger systems, geometry be- 
comes less important and for the same filling we observe 
a spin depolarization as a result of correlations. Finally, 
we compared the fully interacting model with the Hub- 
bard and extended Hubbard models. While qualitative 
agreement for the charge-neutral system was observed, 
the effect of correlations can be described only with the 
inclusion of all direct long-range, exchange, and scatter- 
ing terms. 
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FIG. 13: The excitation gaps corresponding to phase dia- 
grams from Fig. [I2]for many-body Hamiltonians with (a) only 
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support from fellowship cofinanced by European Union 
within European Social Fund. A. W. acknowledges sup- 
port from the EU Marie Curie CIG. 



where 4> z (r) are p z orbitals. The positions of the sub- 
lattice A and B atoms are given by Ra = nai + m&2 
and Rb = na.% + ma.2 + b, described by unit vectors 
of hexagonal lattice defined as ai.2 = a/2(±v / 3, 3) and 
b = a(0, 1), a vector between two nearest-neighbors 
atoms from the same unit cell with a distance a = 1.42 
A. N c is the number of unit cells, and exp(i^k) = |j(k)| 



with /(k) = 1 + e 



ika 2 



The density matrix for the 



graphene layer p° l(7 for two sites j and I is defined as 



k 



(io) 



where 6r's are the coefficients of the p z orbitals given 
in Eq. ©. The on-site density matrix element for an 
arbitrary lattice site j is site and sublattice index inde- 
pendent, 



33° 



2N, 



L E 



e -ikR Je ikR.j 



— E 



1 



1 

2' 



(11) 



where we took into account the fact that the number of 
occupied states is equal to the number of unit cells in the 
system. The nearest-neighbors density matrix elements 
for two atoms from the same unit cell corresponds to 
R; Rj = b and can be calculated using 



Piu, 



— E 

O.N ' 



e -ikRj gikR, g — ikbg- 



0.262, 



where the summation over occupied valence states is car- 
ried out numerically. We obtained the same value for 
two other nearest neighbors. The same results can also 
be obtained by diagonalizing a sufficiently large hexago- 
nal graphene quantum dot and by computing the density 
matrix elements for two nearest neighbors in the vicinity 
of the center of the structure. We have also calculated 
next-nearest neighbors density matrix elements, obtain- 
ing negligibly small values. 



APPENDIX A 

In this appendix, we calculate density matrix elements 
p° l(7 , between sites j and / for an infinite graphene sheet. 
The valence band eigenfunctions of the TB Hamiltonian 
in the nearest-neighbor approximation given by Eq. ([2]) 
are 

n(r) = -^=£> 4kRA <Mr-R A ) 

V c R A 

+ ^ e ^R Be -ikb e -ie k ^ (r _ RB))) (9) 

Rb 
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